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ABSTRACT 


In  recent  years  many  models  of  three-dimensional  (3-D)  seismic  velocity  structure  in  Eurasia  have  been  developed 
using  a  variety  of  techniques  and  data.  Most  of  these  models  are  not  accompanied  by  quantitative  estimates  of 
uncertainty,  either  in  the  model  parameters  themselves  (e.g.  seismic  velocities)  or  in  geophysical  observables 
predicted  by  the  models  (e.g.  body-wave  travel  times).  Moreover,  the  various  3-D  models  produced  by  these  studies 
have  not  been  compared  to  one  another  in  their  predictive  capabilities  in  any  meaningful  way  within  the  monitoring 
research  community.  To  address  these  issues  we  are  developing  and  applying  evaluation  metrics  that  robustly 
quantify  and  compare  the  uncertainty  and  predictive  capability  of  3-D  seismic  velocity  models. 

There  are  two  major  elements  in  our  approach.  First,  we  are  performing  a  comprehensive  evaluation  of  a  set  of  3-D 
velocity  models  for  Eurasia  based  on  previously  developed  data  misfit  and  event  mislocation  metrics.  These  metrics 
are  computed  using  a  ground-truth  (GT)  data  set  available  at  the  International  Seismic  Centre  (ISC).  We  discuss 
some  of  the  progress  we  have  made  in  reconciling  discrepancies  between  the  parameterization  and  coverage  of 
various  models  available  to  the  research  community,  and  compare  some  of  the  standard  metrics  we  have  been  able 
to  derive  for  these  models. 

Second,  we  are  investigating  a  new  approach  to  evaluating  velocity  models  based  on  a  Bayesian  framework  for 
model  uncertainty  in  the  travel-time  tomography  problem.  In  this  approach,  geostatistical  parameters  describing 
velocity  heterogeneity  in  the  Earth  are  estimated  from  travel-time  residuals  observed  along  a  set  of  event-station 
paths.  The  parameters  include  the  velocity  variance,  which  quantifies  the  strength  of  velocity  heterogeneity,  and 
vertical  and  horizontal  correlation  lengths,  which  quantify  the  spatial  smoothness  of  velocity  heterogeneity.  The 
inferred  geostatistical  parameters  characterize  the  discrepancy  between  the  Earth's  velocity  and  the  velocity  of  the 
reference  model  used  to  calculate  the  travel-time  residuals.  That  is,  they  quantify  the  uncertainty  of  the  reference 
model  and  thus  serve  as  metrics  for  model  evaluation.  In  addition,  as  we  have  shown  in  previous  projects, 
geostatistical  parameters  can  be  used  to  calculate  the  uncertainty  in  travel  times  predicted  by  the  reference  model  for 
arbitrary  paths.  We  formulate  the  geostatistical  estimation  problem  within  a  maximum-likelihood  framework  and 
outline  an  approach  for  its  numerical  solution. 

We  are  currently  pursuing  a  maximum  likelihood  approach  to  solving  the  geostatistical  estimation  problem  whereby 
the  geostatistical  parameters,  together  with  pick-error  variances,  are  simultaneously  fit  to  the  second-order  statistics 
of  the  observed  travel-time  residuals.  The  maximum  likelihood  criterion  for  the  optimal  parameter  values  reduces  to 
a  set  of  coupled,  nonlinear  equations  which,  in  general,  must  be  solved  numerically.  We  formulate  the  geostatistical 
estimation  problem  within  a  maximum-likelihood  framework  and  outline  an  approach  for  its  numerical  solution. 
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OBJECTIVES 

The  main  objective  of  our  project  is  to  develop  meaningful  measures  of  the  predictive  capabilities  of  the  multitude 
of  models  that  have  been  developed  over  recent  years  by  the  nuclear  monitoring  research  community.  While  our 
primary  focus  is  on  regional-scale  seismic  velocity  models,  we  intend  our  general  approach  to  be  applicable  as  well 
to  global  Earth  models  and  models  of  other  geophysical  parameters  such  as  attenuation  and  density. 

Our  project  consists  of  two  major  elements.  First,  we  have  collected  a  set  of  regional  3-D  velocity  models  for  Asia 
and  are  performing  a  comprehensive  and  methodical  evaluation  of  them  using  data  misfit  and  event  mislocation 
metrics.  Metrics  will  be  evaluated  for  each  model  based  on  a  common  ground-truth  data  set  comprising  GT5  events 
and  their  arrival  picks.  A  large  part  of  this  effort  so  far  has  been  spent  developing  standardized  methods  for  model 
parameterization,  forward  modeling  (e.g.  ray  tracing),  and  event  relocation.  Some  of  the  tools  we  need  for  these  tests 
are  available  to  us  in-house  from  our  previous  tomography  and  location  projects,  but  others  are  being  gathered  from 
outside  sources. 

The  second  element  of  our  project  is  the  investigation  of  a  new  approach  for  estimating  the  uncertainty  and 
predictive  capability  of  3D  velocity  models.  Our  approach  will  quantify  model  uncertainty  in  terms  of  spatially 
variable  geostatistical  parameters  fit  to  second-order  statistics  of  travel-time  residuals,  and  then  convert  these 
parameters  into  the  uncertainty  in  travel-time  or  other  geophysical  predictions.  To  do  this  requires  covariance 
modeling  capabilities  we  developed  for  travel  times  in  an  earlier  project  (Rodi  and  Myers,  2008)  and  which  can  be 
extended  to  other  observables. 

RESEARCH  ACCOMPLISHED 


3-D  Seismic  Velocity  Models,  Ground-Truth  Data,  and  Analysis  Techniques  to  Examine  Model  Performance 

Our  first  task  has  been  to  collect  3-D  regionalized  seismic  velocity  models  and  format  them  in  a  way  that  makes 
comparisons  straightforward.  Table  1  lists  the  six  velocity  models  that  we  are  currently  using  in  our  study,  along 
with  their  respective  authors  and  a  brief  description.  Our  test-bed  model  for  the  project  is  the  Joint  Weston/MIT 
(JWM)  inversion  model  (Reiter  and  Rodi,  2009).  JWM  is  the  result  of  a  joint  inversion  of  regional  P  travel  times 
and  Rayleigh  fundamental-mode  group  velocities.  It  consists  of  the  P  and  S  velocities  and  density  of  the  crust  and 
upper  mantle  for  a  broad  region  of  Asia  (defined  in  a  geographic  box  from  10  -  50°N,  40  -  1 10°E).  In  addition  to 
the  JWM  model,  we  have  obtained  the  CUB2.0  surface-wave  inversion  model  (Ritzwoller  et  al.,  2002)  with  density; 
the  Stevens  et  al.  (2008)  surface-wave  inversion  model;  the  EAV09  inversion  model  from  the  group  headed  by 
Suzan  van  der  Lee  at  Northwestern  University  (Schmid  et  al.,  2008);  and  a  new  P-velocity  inversion  model 
(LLNLG3-D)  from  Nathan  Simmons  and  Stephen  Myers  (Simmons  et  al.,  2011).  We  also  have  access  to  an  a 
priori  regionalized  model  known  as  the  DOE  Unified  Model  (Begnaud  et  al.,  2004;  Pasyanos  et  al.,  2004).  We 
continue  to  seek  models  from  other  researchers,  but  many  are  not  openly  available  or  do  not  meet  the  necessary 
criteria  to  be  included  in  our  study. 

One  of  the  more  difficult  problems  we  have  encountered  in  the  project  has  been  the  difficulty  in  accurately 
representing  the  models  in  our  study  using  a  single  common  format.  At  the  present  time  there  is  no  standardized 
format  for  model  exchange,  which  has  prevented  the  community  and  program  funding  managers  from  understanding 
the  significant  differences  between  3-D  models  that  are  constructed  for  use  in  predicting  nuclear  monitoring 
observables.  Our  in-house  model  format  represents  a  3-D  model  as  a  geographic  grid  of  1-D  profiles,  with  each 
profile  sharing  a  common  parameterization  with  respect  to  depth.  For  example,  in  the  JWM  model  the  crust  at  each 
latitude/longitude  point  is  divided  into  vertically  homogeneous  layers  corresponding  to  the  those  of  the  CRUST2.0 
model:  water  overlying  two  sediment  layers  (soft  and  hard  sediments)  overlying  three  metamorphic/igneous  layers 
(upper,  middle  and  lower  crust).  The  upper  mantle  is  represented  vertically  by  piecewise  linear  parameter  functions 
sampled  at  multiple  nodes  distributed  between  the  Moho  discontinuity  and  a  depth  of  410  km.  The  versatility  of  this 
model  format  allows  discontinuities  and  homogeneous  layers  to  be  accurately  captured,  and  no  special  software  is 
necessary  to  allow  the  models  to  be  exchanged  between  researchers.  Most  of  the  models  in  our  study  can  be 
represented  in  this  format,  although  some  conversion  is  necessary  in  certain  cases  (LP2008;  Stevens  et  al.,  2008). 

We  believe  that  the  seismic  monitoring  community  would  be  well  served  by  developing  and  adopting  a  standard 
Earth  model  exchange  format.  Other  research  communities  have  long-standing  efforts  in  this  area.  For  example,  the 
international  climate  modeling  community  developed  a  standard  protocol  in  1990  for  global  atmospheric  general 
circulation  models  (AGCMs)  that  provides  a  framework  in  support  of  climate  model  diagnosis,  validation. 
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intercomparison,  documentation  and  data  access.  Their  framework  enables  a  diverse  community  of  scientists  to 
analyze  AGCMs  in  a  systematic  fashion,  which  serves  to  facilitate  model  improvement  (Gates,  1992).  Our  research 
community  would  also  benefit  greatly  from  the  development  of  a  framework  for  model  exchange,  similar  to  the  one 
the  seismic  community  has  for  data  exchange. 

Table  1.  Descriptions  of  the  current  set  of  3-D  seismic  velocity  models  used  in  the  study. 


Model  Name 

Authors 

Geographic 

Coverage 

Type 

References 

JWM 

Weston/MIT 

10  -  50°N, 

30  -  120°  E 

inversion  model  of  Vp,  Vs  and  density; 
constructed  with  Pn  travel  times  and  Rayleigh 
group  velocities 

Reiter  and  Rodi  (2009) 

LLNLG3D 

LLNL 

Global,  high- 

res  over 

MidEast 

inversion  model  of  Vp  constructed  fromP  and  Pn 
travel  times 

Myers  et  al.  (201 1); 
Simmons  et  al .  (201 1) 

DOE  Unified  Model 

LLNI7LANL 

0  -  85°N,  20 

-  75°  E 

a  priori  model  of  Vp,  Vs,  density,  Qp  and  Qs 

Pasyanos  et  al .,  2004; 
Begnaud  et  al .,  2004 

CUB20_J362D28 

CU  Boulder 

Harvard 

Global 

inversion  model  of  Vp,  Vs  and  density  fromO  -  700 
km  defined  over  the  globe;  determined  with  group 
velocities 

Shapiro  and  Ritzwoller, 
2002 

LP2008 

SAIC 

Global 

inversion  model  of  Vp,  Vs,  density,  Qbeta  fromO  - 
—700  km  defined  over  the  globe;  determined  with 
surface  waves 

Stevens  et  al .,  2008 

EAV09 

Northwestern 

Univ./  LLNL 

10  -  60°N, 

-35  -  80°  E 

inversion  model  of  Vs  with  add'l  conversion  to  Vp; 
determined  mainly  from  waveforms  and  teleseismic 
S  travel  times 

Schmid  et  al .,  2008 

To  demonstrate  the  usefulness  of  performing  simple  model  comparisons,  Figure  1  displays  map-view  slices  from  the 
models  in  our  study.  The  geographic  coverage  of  the  maps  in  Figure  1  was  chosen  based  on  the  footprint  of  the 
JWM  model,  which  means  that  two  of  the  models  do  not  completely  cover  their  respective  plots  (LLNL  G3D  is 
global,  but  we  only  requested  a  certain  piece  of  it,  and  the  EAV09  model  is  concentrated  on  the  Middle  East).  The 
Generic  Mapping  Tool  (GMT)  package  (Wessel  and  Smith,  1995)  was  used  to  plot  the  P  and  S  velocities  at  a  depth 
of  90  km,  using  the  same  color  map  for  each  model.  Prior  to  creating  this  plot,  the  models  were  reformatted  to  our 
in-house  3-D  parameterization  and  then  interpolated  to  a  depth-slice  grid.  The  different  panels  in  Figure  1  offer  a 
direct  comparison  of  the  models  and  illustrate  their  similarities  and  differences.  For  example,  JWM,  LLNL  G3D 
and  CUB2.0  are  similar  to  each  other  in  the  strength  of  heterogeneities  at  this  upper-mantle  depth.  All  models 
demonstrate  similar  patterns  of  slow-  and  fast-velocity  areas  across  the  plotted  region.  CUB2.0  is  defined  on  a  2°x2° 
grid,  which  produces  a  smoother  result  compared  to  the  other  models.  The  EAV09  model  is  similar  to  other  models 
in  its  shear  velocity  values,  but  the  P  velocities  appear  to  exhibit  weaker  overall  variation  compared  to  the  other 
models.  The  LP2008  model  is  ‘blocky’  in  nature,  which  is  to  be  expected  since  it  is  comprised  of  homogeneous 
layers  that  are  conducive  to  dispersion  modeling  calculations.  The  a  priori  LLNL  Unified  model  is  least  like  the 
others  in  pattern  and  heterogeneity  strength  -  we  note  that  this  model  is  equivalent  to  the  3SMAC  model  (Nataf  and 
Ricard,  1996)  at  this  depth. 
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Depth  Slices  at  Z  =  90  km 


Vp  model  only: 
Vs  not  defined 


7.6$  7  JO  7  .as  ,410  6.25  a  JO 
Vp(krrVs) 


7.65  7.80  7.95  8.10  825  840 
Vp  (km/s) 

*  5<f  60’  70°  aor  90’  iocr  nor 


4,24  4.32  4.40  448  4.56  464 
Vs(km/s) 


LLNLG3D 


7.65  7 50  755  8.10  825  840 
Vp  (km/s) 


424  422  4.40  448  456  464 
Vs  (km/s) 


7.65  750  755  8.10  825  840 
Vp  (km/s) 


424  422  440  448  456  464 
Vs  (km/s) 


CUB2 


EAV09 


7.65  7.80  755  8.10  825  840 
Vp  (km/s) 


DOE  Unified  Model 


4.24  422  440  448  456  464 
Vs  (km/s) 


LP2008 


Figure  1.  Map-view  slices  of  P  and  S  velocity  at  90  km  depth  from  the  regional  3-D  velocity  models  in  our 

study  (see  Table  1).  The  models  are  labeled  at  the  bottom  of  each  Vp  and  Vs  plot,  and  the  color  scale 
is  the  same  across  all  panels.  See  the  text  for  further  explanation  of  these  plots. 

To  evaluate  models  with  respect  to  one  another,  we  require  a  high-quality  set  of  GT  events  whose  epicentral 
locations  are  known  precisely.  We  used  the  new  catalog  of  GTO-5  reference  events  (Bondar  and  McLaughlin,  2009) 
hosted  by  the  ISC  (http://www.isc.ac.uk).  This  global  database  includes  more  than  7,000  events  whose  epicentral 
location  accuracy  is  known  to  at  least  5  km.  GT  events  with  well-established  locations  and  origin  times  have  been 
used  by  multiple  authors  to  validate  velocity  models  (Ritzwoller  et  al.,  2003;  Flanagan  et  al.,  2007).  We  filtered  the 
new  International  Association  of  Seismology  and  Physics  of  the  Earth's  Interior  (IASPEI)  Reference  Event  List 
(REL)  for  events  in  our  study  region  and  found  248  GTO-5  explosions  and  348  GT5  earthquakes,  most  of  which  are 
in  specific  geographic  clusters.  From  these  events  we  extracted  the  defining  regional  P  and  S  arrival-time  picks  in 
the  IASPEI  REL  that  were  used  to  develop  the  GT  locations. 

Using  a  variety  of  filtering  criteria  designed  to  eliminate  outliers,  we  derived  a  validation  data  set  of  9,242  P-wave 
and  2,214  S-wave  regional  arrivals  observed  at  stations  within  a  latitude-longitude  box  defined  by  (0-60°  N,  30-120° 
E).  The  great-circle  raypath  coverage  for  the  P  and  S  GT  observations  is  shown  in  Figure  2.  While  the  S  GT 
observations  are  not  strictly  needed  for  most  of  the  model  validation  exercises  that  we  perform,  we  included  them 
with  the  idea  they  may  be  useful  in  future  studies.  It  is  clear  that  the  current  version  of  our  GT  data  ray  paths  sample 
only  a  limited  portion  of  the  JWM  study  region,  which  illustrates  the  difficulty  of  validating  a  model  with  travel 
times  alone.  However,  the  IASPEI  REL  database  is  currently  the  highest  quality  GT  data  resource  we  have  to 
demonstrate  the  performance  of  the  3-D  models. 
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▲  Stations  o  Reference  Events 

Figure  2.  Great-circle  raypath  coverage  for  the  P  (left)  and  S  (right)  paths  in  the  IASPEI REL  ground-truth 
database.  Stations  are  represented  by  purple  triangles  and  events  by  gray  circles.  Note  the  sparse 
coverage  over  the  area  of  the  regional  3-D  models. 

We  have  begun  to  augment  the  GTO-5  IASPEI  REL  data  set  with  other  observations  that  have  been  collected  in 
previous  research  efforts.  While  these  supplemental  databases  may  only  be  of  GT 10-20  quality  in  the  event 
locations,  they  can  still  provide  insights  into  the  travel-time  prediction  behavior  of  the  various  models,  and  they  are 
a  valuable  data  source  for  testing  the  variance  estimation  approach  we  are  pursuing  under  the  second  major  task  of 
the  project.  For  instance,  a  consortium  effort  led  by  SAIC  in  the  early  2000s  produced  several  valuable  GT  data 
resources  that  are  not  included  in  the  IASPEL  REL  (Murphy  et  ah,  2005).  The  electronic  supplement  that 
accompanies  the  Murphy  et  al.  (2005)  paper  contains  several  groomed  event  bulletins  that  we  plan  to  utilize  in  our 
study.  One  of  these  is  specifically  focused  on  China  and  consists  of  travel-time  observations  at  49  stations  of  the 
dense  Chinese  National  Network  that  were  used  in  a  tomographic  inversion.  Murphy  et  al.  (2005)  notes  that  the 
average  quality  of  the  Chinese  bulletin  locations  is  on  the  order  of  GT10,  making  them  suitable  for  use  in  a  broad 
number  of  validation  studies.  Figure  3  shows  the  distribution  of  stations  and  events  in  the  Chinese  tomographic 
inversion  bulletin  that  we  have  begun  to  analyze  for  use  in  metric  performance  calculations. 
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Figure  3.  Supplemental  ground-truth  data  set  from  the  Chinese  ABCE  bulletin,  collected  under  a  previous 

nuclear  monitoring  research  consortium  effort  (Murphy  et  al.,  2005).  The  events  are  shown  as  black 
dots,  and  the  Chinese  stations  as  purple  triangles. 
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While  Figure  1  illustrates  that  the  3-D  models  in  our  study  can  seem  quite  different  in  snapshot  form,  they  often 
perform  like  each  other  in  travel-time  prediction  and  location  tests.  To  examine  these  behaviors  we  are  investigating 
the  performance  of  the  3-D  models  using  different  forward  modeling  techniques,  in  an  attempt  to  determine 
whether,  how  and  why  they  differ.  For  example,  in  Figure  4  we  show  the  results  of  predicting  the  ray  path  from  a 
GTO  explosion  in  the  Western  Soviet  Union  to  station  SHI  in  Iran,  using  all  six  models  in  Table  1  and  the  Podvin- 
Lecomte  (1991)  ray  tracing  method.  The  epicentral  path  length  for  the  event-station  pair  is  approximately  13.6°  and 
in  the  IASPEL  REL  the  bulletin  observation  is  194  s  for  the  travel  time  of  the  Pn  phase.  The  station-path  map  is 
shown  in  the  center  of  Figure  4,  and  arrayed  around  the  map  are  summed  projections  (collapsed  to  Cartesian 
latitude)  of  the  ray  paths  that  each  3-D  model  produces  for  the  station-event  pair  that  we  chose.  In  each  inset  ray 
path  sub-panel  we  list  the  bottoming  point  of  the  ray  and  the  travel  time  along  the  path. 


The  results  show  that  three  of  the  models  produce  ray  paths  that  meet  a  classical  definition  of  the  Pn  phase,  traveling 
near  the  base  of  the  Moho  discontinuity  in  the  uppermost  mantle.  However,  three  of  the  models  produce  ray  paths 
that  dive  deeper  into  the  mantle  along  a  traditional  P-wave  trajectory,  bottoming  at  depths  greater  than  140  km.  The 
path  travel  times  also  vary  fairly  significantly  across  the  six  models,  with  a  spread  of  ~7  seconds.  It  is  interesting  to 
note  that  all  of  the  models  have  Moho  depths  close  to  45  km  on  both  the  station  and  event  sides  of  the  ray  paths,  so 
further  analysis  is  warranted  to  explain  the  wide  variation  in  the  ray  path  behaviors  for  this  case. 
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Figure  4.  Example  illustrating  the  ray  path  behavior  in  each  3-D  model  for  a  single  event-to-station  pair  in 
the  ground-truth  database.  The  rays  are  comprised  of  the  sensitivities  calculated  along  the 
particular  path  by  the  Podvin-Lecomte  method;  they  are  depicted  as  projections  (summations)  to 
latitude. 


The  results  of  this  type  of  exercise  illustrate  the  significant  differences  that  can  exist  in  the  subsurface  spatial 
sampling  and  travel  times  of  a  given  ray  path  across  the  various  models.  We  are  analyzing  the  differences  between 
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different  ray  tracing  techniques  (e.g.  ray-bending,  ray  shooting,  fast  marching  method,  eikonal)  to  see  how  the 
model  performance  changes.  Our  aim  is  to  determine  whether  the  prediction/ forward  modeling  techniques  used  to 
develop  individual  3-D  Earth  models  must  be  included  as  part  of  the  distribution  of  a  new  model  within  the  research 
community. 


Evaluation  of  Velocity  Models  Based  on  Model  Uncertainty  Analysis 

To  address  the  second  major  task  in  our  project,  we  have  developed  a  mathematical  framework  for  the  problem  of 
estimating  geostatistical  parameters  of  velocity  heterogeneity  in  the  Earth  from  travel-time  residuals  observed  along 
a  set  of  event-station  paths.  The  inferred  geostatistical  parameters  characterize  the  discrepancy  between  the  Earth's 
velocity  and  the  velocity  of  the  reference  model  used  to  calculate  the  travel-time  residuals.  That  is,  they  quantify  the 
uncertainty  of  the  reference  model  and  thus  serve  as  metrics  for  model  validation.  In  addition,  as  we  have  shown  in 
previous  projects,  geostatistical  parameters  can  be  used  to  calculate  the  uncertainty  in  travel  times  predicted  by  the 
reference  model  for  arbitrary  paths. 

The  problem  can  be  formulated  mathematically  as  follows.  Let  the  vector  m  contain  parameters  describing  the 
difference  between  the  slowness  functions  of  the  real  Earth  and  a  reference  model.  Then  a  vector  of  observed  travel¬ 
time  residuals  r,  calculated  with  respect  to  the  reference  model,  can  be  related  to  m  to  first  order  by 

r  =  Am  +  e,  (1) 

where  A  is  a  travel-time  sensitivity  matrix  and  e  is  a  vector  containing  measurement  (pick)  errors  in  the  residuals. 
The  pick  errors  are  assumed  to  be  zero  mean  with  some  variance/covariance  matrix  Ce.  Likewise,  we  assume  m  is 
zero  mean  with  variance  C„,.  Equation  1  then  implies  that  r  is  zero  mean  with  variance  matrix  given  by 

E[rr7>  Cr  =  ACmAT  +  Ce.  (2) 

Now  assume  that  Ce  is  a  given  function  of  a  parameter  vector  0  =  (0X  02  ■■■Y ,  and  that  C,„  is  parameterized  by 
Y  —  (Yi  Vi  "•  Y  ■  The  0 for  example,  might  be  pick-error  variances  as  a  function  of  epicentral  distance.  The  yk  are 
the  geostatistical  parameters  of  interest,  comprising  slowness  variances  and  correlation  lengths  as  a  function  of 
position  in  the  Earth.  The  problem  at  hand  is  to  infer  the  parameters  9k  and  yk  from  r  on  the  basis  of  Equation  1, 
recognizing  that  the  yk  are  the  primary  targets. 


We  have  examined  a  number  of  possible  approaches  to  solving  this  variance  estimation  problem  and  are  currently 
pursuing  a  formal  approach  based  on  maximum-likelihood  estimation.  The  relevant  likelihood  function  is  provided 
by  the  probability  density  function  of  r.  Assuming  m  and  e,  and  thus  r,  are  Gaussian,  the  likelihood  function,  L,  is 
given  by 

logL  ^0,  y)  =  const  —  i[log(detCr)  +  rrCiT1r],  (3) 


where  it  is  understood  that  Cr  depends  on  8  and  y.  The  optimal  estimates  of  these  parameter  vectors  are  the  values 
that  maximize  L.  This  implies  that  L  is  stationary  with  respect  to  the  parameters,  or 


(4) 


Given  Equations  2  and  3,  and  the  properties  of  matrix  determinants,  the  stationarity  conditions  become 

r  (5) 

(6) 


trace  Cr  1  —  =  rrCr  1  —  Cr  1 


trace  C^A  —  Ar  =  r^Cj^A— ATC“1r. 

r  aYk  r  dYk  1 

These  equations  can  be  simplified  by  substituting  the  theoretical  results  of  tomographic  inversion  implied  by  the 
stochastic  assumptions  we  have  made.  That  is,  the  maximum  a  posteriori  (MAP)  estimate  of  m,  and  its  residual 
vector,  are  given  by 


m  =  CmA‘  Cr1r 


e  =  r  -  Am  =  CeCr  1r. 


(7) 

(8) 
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Also  useful  are  the  influence  and  resolution  matrices,  S  and  R,  respectively: 


S  =  AC^C’1  (9) 

R  =  CmATCr 1A.  (10) 

Equations  5  and  6  can  then  be  rewritten  as 

trace  (I-Sjgcj1  =  e^CJ1  (11) 

trace  R^C”1  =  m7'C“1  ^C^m.  (12) 


Equations  11  and  12  are  coupled,  nonlinear  equations  that  cannot  be  solved  analytically  except  in  very  simple 
situations.  In  a  problem  involving  a  large  number  of  residuals  and  slowness  parameters,  efficient  numerical  solution 
schemes  are  also  elusive.  The  difficulty  of  the  problem  depends  on  the  particular  choice  of  the  parameters  0 \ and 
We  are  in  the  process  of  devising  iterative,  numerical  algorithms  for  solving  Equations  1 1  and  12  for  pick-error  and 
slowness  variances,  holding  correlation  lengths  of  slowness  heterogeneity  fixed.  If  the  resulting  algorithms  are 
efficient,  the  solution  for  correlation  lengths  can  possibly  be  found  by  grid  search  or  trial-and-error,  if  more  direct 
schemes  are  not  discovered. 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  main  objective  of  our  research  project  is  to  develop  and  apply  meaningful  measures  of  the  predictive 
capabilities  of  3-D  Earth  models.  For  our  first  major  task,  we  are  performing  a  comprehensive  and  methodical 
evaluation  of  a  set  of  regional  velocity  models  for  Asia  based  on  data  misfit  and  event  mislocation  metrics.  We  are 
currently  working  with  six  3-D  regional  velocity  models  of  the  crust  and  mantle  and  are  specifically  developing  our 
techniques  so  that  they  are  easily  applicable  to  other  models  that  may  become  available.  In  the  past  year  our  work  on 
this  first  task  has  focused  on  the  accurate  conversion  of  models  from  their  native  formats  into  a  single  in-house 
representation  and  the  testing  of  the  numerical  tools  (such  as  a  set  of  ray  tracing  methods)  to  predict  travel  times 
using  the  models. 

In  our  second  task  we  are  developing  a  new  method  to  evaluate  models  based  on  the  analysis  of  their  uncertainty. 

We  are  currently  pursuing  a  Bayesian  approach  to  this  problem  in  which  geostatistical  parameters  describing 
velocity  uncertainty,  together  with  pick-error  variances,  are  simultaneously  fit  to  the  second-order  statistics  of 
observed  travel-time  residuals.  To  date,  we  have  formulated  theoretical  solutions  within  the  framework  of 
maximum-likelihood  estimation  and  are  seeking  practical  numerical  techniques  for  computing  such  solutions. 
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